﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Runtime.Serialization;
using Jama.Common;
using Jama.SVD;
using Jama.LU;
using Jama.QR;
using Jama.Cholesky;
using Jama.Eigenvalue;

namespace Jama
{
    public class Matrix //: ICloneable, ISerializable
    {
        #region Attributes
        /** Array for internal storage of elements.
        @serial internal array storage.
        */
        private double[][] A;

        /** Row and column dimensions.
        @serial row dimension.
        @serial column dimension.
        */
        private int m, n;
        #endregion

        #region Constructor(s)
        /** Construct an m-by-n matrix of zeros. 
        @param m    Number of rows.
        @param n    Number of colums.
        */
        public Matrix (int m, int n)
        {
          this.m = m;
          this.n = n;
          A = new double[m][];
          
            for(int i = 0 ; i < m ; i++ )    
              A[i] = new double[n];
        }

        /** Construct an m-by-n constant matrix.
        @param m    Number of rows.
        @param n    Number of colums.
        @param s    Fill the matrix with this scalar value.
        */
        public Matrix (int m, int n, double s) 
        {
            this.m = m;
            this.n = n;
            A = new double[m][];
            for (int i = 0; i < m; i++) 
            {
                A[i]= new double[n];
                for (int j = 0; j < n; j++)
                {
                A[i][j] = s;
                }
            }
        }

        /** Construct a matrix from a 2-D array.
        @param A    Two-dimensional array of doubles.
        @exception  IllegalArgumentException All rows must have the same length
        @see        #constructWithCopy
        */
        public Matrix (double[][] A) 
        {
            m = A.Length;
            n = A[0].Length;
            for (int i = 0; i < m; i++)
            {
                if (A[i].Length != n) 
                {
                throw new IllegalArgumentException("All rows must have the same length.");
                }
            }
            this.A = A;
        }

        /** Construct a matrix quickly without checking arguments.
        @param A    Two-dimensional array of doubles.
        @param m    Number of rows.
        @param n    Number of colums.
        */
        public Matrix (double[][] A, int m, int n)
        {
            this.A = A;
            this.m = m;
            this.n = n;
        }

        /** Construct a matrix from a one-dimensional packed array
        @param vals One-dimensional array of doubles, packed by columns (ala Fortran).
        @param m    Number of rows.
        @exception  IllegalArgumentException Array length must be a multiple of m.
        */
        public Matrix (double[] vals, int m)
        {
            this.m = m;
            n = (m != 0 ? vals.Length/m : 0);
            if (m*n != vals.Length) 
            {
                throw new IllegalArgumentException("Array length must be a multiple of m.");
            }
            A = new double[m][];
            for (int i = 0; i < m; i++)
            {
                A[i] = new double[n];
                for (int j = 0; j < n; j++)
                {
                A[i][j] = vals[i+j*m];
                }
            }
        }
        #endregion

        #region Public Methods
        /// <summary>
        /// Construct a matrix from a copy of a 2-D array.
        /// throws IllegalArgumentException All rows must have the same length
        /// </summary>
        /// <param name="A">Two-dimensional array of doubles.</param>
        /// <returns></returns>
        public static Matrix constructWithCopy(double[][] A) 
        {
            int m = A.Length;
            int n = A[0].Length;
            Matrix X = new Matrix(m,n);
            double[][] C = X.getArray();

            for (int i = 0; i < m; i++) 
            {
                if (A[i].Length != n)
                {
                    throw new IllegalArgumentException ("All rows must have the same length.");
                }
                for (int j = 0; j < n; j++)
                {
                    C[i][j] = A[i][j];
                }
            }
            return X;
        }

        /// <summary>
        /// Make a deep copy of a matrix
        /// </summary>
        /// <returns></returns>
        public Matrix copy ()
        {
            Matrix X = new Matrix(m,n);
            double[][] C = X.getArray();
            for (int i = 0; i < m; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    C[i][j] = A[i][j];
                }
            }
            return X;
        }

        /// <summary>
        /// Clone the Matrix object.
        /// </summary>
        /// <returns></returns>
        public Object clone ()
        {
            return this.copy();
        }
        
        /// <summary>
        ///  Access the internal two-dimensional array.   
        /// </summary>
        /// <returns>
        /// Pointer to the two-dimensional array of matrix elements.
        /// </returns>
        public double[][] getArray () 
        {
            return A;
        }

        /// <summary>
        /// Copy the internal two-dimensional array.
        /// </summary>
        /// <returns>
        /// Two-dimensional array copy of matrix elements.
        /// </returns>
        public double[][] getArrayCopy ()
        {
            double[][] C = new double[m][];
            for (int i = 0; i < m; i++)
            {
                C[i]= new double[n];
                for (int j = 0; j < n; j++)
                {
                C[i][j] = A[i][j];
                }
            }
            return C;
        }

        /// <summary>
        ///  Make a one-dimensional column packed copy of the internal array.
        /// </summary>
        /// <returns>
        /// Matrix elements packed in a one-dimensional array by columns.
        /// </returns>
        public double[] getColumnPackedCopy () 
        {
            double[] vals = new double[m*n];
            for (int i = 0; i < m; i++)
            {
                for (int j = 0; j < n; j++)
                {
                vals[i+j*m] = A[i][j];
                }
            }
            return vals;
        }

        /// <summary>
        /// Make a one-dimensional row packed copy of the internal array.
        /// </summary>
        /// <returns>
        /// Matrix elements packed in a one-dimensional array by rows.
        /// </returns>
        public double[] getRowPackedCopy () 
        {
            double[] vals = new double[m*n];
            for (int i = 0; i < m; i++) 
            {
                for (int j = 0; j < n; j++)
                {
                    vals[i*n+j] = A[i][j];
                }
            }
            return vals;
        }
        
        /// <summary>
        /// Get row dimension.
        /// </summary>
        /// <returns>
        /// m, the number of rows.
        /// </returns>
        public int getRowDimension () 
        {
            return m;
        }

        /// <summary>
        /// Get column dimension.
        /// </summary>
        /// <returns>
        /// n, the number of columns.
        /// </returns>
        public int getColumnDimension ()
        {
            return n;
        }

        /// <summary>
        /// Get a single element.
        /// throws ArrayIndexOutOfBoundsException
        /// </summary>
        /// <param name="i">Row index.</param>
        /// <param name="j"> Column index.</param>
        /// <returns>A(i,j)</returns>
        public double get (int i, int j)
        {
            return A[i][j];
        }

   
        /// <summary>
        /// Get a submatrix.
        /// throws ArrayIndexOutOfBoundsException Submatrix indices
        /// </summary>
        /// <param name="i0">Initial row index</param>
        /// <param name="i1">Final row index</param>
        /// <param name="j0">Initial column index</param>
        /// <param name="j1">Final column index</param>
        /// <returns>A(i0:i1,j0:j1)</returns>
        public Matrix getMatrix (int i0, int i1, int j0, int j1) 
        {
            Matrix X = new Matrix(i1-i0+1,j1-j0+1);
            double[][] B = X.getArray();
            try 
            {
                for (int i = i0; i <= i1; i++) 
                {
                for (int j = j0; j <= j1; j++)
                {
                    B[i-i0][j-j0] = A[i][j];
                }
                }
            } catch(Exception ex)
            {
                throw new ArrayIndexOutOfBoundsException("Submatrix indices");
            }
            return X;
        }

   /** Get a submatrix.
   @param r    Array of row indices.
   @param c    Array of column indices.
   @return     A(r(:),c(:))
   @exception  ArrayIndexOutOfBoundsException Submatrix indices
   */

   public Matrix getMatrix (int[] r, int[] c)
   {
      Matrix X = new Matrix(r.Length,c.Length);
      double[][] B = X.getArray();
      try 
      {
         for (int i = 0; i < r.Length; i++) 
         {
            for (int j = 0; j < c.Length; j++)
            {
               B[i][j] = A[r[i]][c[j]];
            }
         }
      } catch(Exception e) 
      {
         throw new ArrayIndexOutOfBoundsException("Submatrix indices");
      }
      return X;
   }

   /** Get a submatrix.
   @param i0   Initial row index
   @param i1   Final row index
   @param c    Array of column indices.
   @return     A(i0:i1,c(:))
   @exception  ArrayIndexOutOfBoundsException Submatrix indices
   */

   public Matrix getMatrix (int i0, int i1, int[] c)
   {
      Matrix X = new Matrix(i1-i0+1,c.Length);
      double[][] B = X.getArray();
      try 
      {
         for (int i = i0; i <= i1; i++)
         {
            for (int j = 0; j < c.Length; j++)
            {
               B[i-i0][j] = A[i][c[j]];
            }
         }
      } catch(Exception e) 
      {
         throw new ArrayIndexOutOfBoundsException("Submatrix indices");
      }
      return X;
   }

   /** Get a submatrix.
   @param r    Array of row indices.
   @param i0   Initial column index
   @param i1   Final column index
   @return     A(r(:),j0:j1)
   @exception  ArrayIndexOutOfBoundsException Submatrix indices
   */

   public Matrix getMatrix (int[] r, int j0, int j1) 
   {
      Matrix X = new Matrix(r.Length,j1-j0+1);
      double[][] B = X.getArray();
      try 
      {
         for (int i = 0; i < r.Length; i++)
         {
            for (int j = j0; j <= j1; j++)
            {
               B[i][j-j0] = A[r[i]][j];
            }
         }
      } catch(Exception e)
      {
         throw new ArrayIndexOutOfBoundsException("Submatrix indices");
      }
      return X;
   }

   /** Set a single element.
   @param i    Row index.
   @param j    Column index.
   @param s    A(i,j).
   @exception  ArrayIndexOutOfBoundsException
   */

   public void set (int i, int j, double s)
   {
      A[i][j] = s;
   }

   /** Set a submatrix.
   @param i0   Initial row index
   @param i1   Final row index
   @param j0   Initial column index
   @param j1   Final column index
   @param X    A(i0:i1,j0:j1)
   @exception  ArrayIndexOutOfBoundsException Submatrix indices
   */

   public void setMatrix (int i0, int i1, int j0, int j1, Matrix X) 
   {
      try
      {
         for (int i = i0; i <= i1; i++)
         {
            for (int j = j0; j <= j1; j++) 
            {
               A[i][j] = X.get(i-i0,j-j0);
            }
         }
      } catch(Exception e)
      {
         throw new ArrayIndexOutOfBoundsException("Submatrix indices");
      }
   }

   /** Set a submatrix.
   @param r    Array of row indices.
   @param c    Array of column indices.
   @param X    A(r(:),c(:))
   @exception  ArrayIndexOutOfBoundsException Submatrix indices
   */

   public void setMatrix (int[] r, int[] c, Matrix X) 
   {
      try
      {
         for (int i = 0; i < r.Length; i++) 
         {
            for (int j = 0; j < c.Length; j++)
            {
               A[r[i]][c[j]] = X.get(i,j);
            }
         }
      } catch(Exception e) 
      {
         throw new ArrayIndexOutOfBoundsException("Submatrix indices");
      }
   }

   /** Set a submatrix.
   @param r    Array of row indices.
   @param j0   Initial column index
   @param j1   Final column index
   @param X    A(r(:),j0:j1)
   @exception  ArrayIndexOutOfBoundsException Submatrix indices
   */

   public void setMatrix (int[] r, int j0, int j1, Matrix X)
   {
      try
      {
         for (int i = 0; i < r.Length; i++) 
         {
            for (int j = j0; j <= j1; j++)
            {
               A[r[i]][j] = X.get(i,j-j0);
            }
         }
      } catch(Exception e)
      {
         throw new ArrayIndexOutOfBoundsException("Submatrix indices");
      }
   }

   /** Set a submatrix.
   @param i0   Initial row index
   @param i1   Final row index
   @param c    Array of column indices.
   @param X    A(i0:i1,c(:))
   @exception  ArrayIndexOutOfBoundsException Submatrix indices
   */

   public void setMatrix (int i0, int i1, int[] c, Matrix X)
   {
      try 
      {
         for (int i = i0; i <= i1; i++)
         {
            for (int j = 0; j < c.Length; j++) 
            {
               A[i][c[j]] = X.get(i-i0,j);
            }
         }
      } catch(Exception e)
      {
         throw new ArrayIndexOutOfBoundsException("Submatrix indices");
      }
   }

   /** Matrix transpose.
   @return    A'
   */

   public Matrix transpose () 
   {
      Matrix X = new Matrix(n,m);
      double[][] C = X.getArray();
      for (int i = 0; i < m; i++) 
      {
         for (int j = 0; j < n; j++)
         {
            C[j][i] = A[i][j];
         }
      }
      return X;
   }

   /** One norm
   @return    maximum column sum.
   */

   public double norm1 () 
   {
      double f = 0;
      for (int j = 0; j < n; j++)
      {
         double s = 0;
         for (int i = 0; i < m; i++) 
         {
            s += Math.Abs(A[i][j]);
         }
         f = Math.Max(f,s);
      }
      return f;
   }

   /** Two norm
   @return    maximum singular value.
   */

   public double norm2 () 
   {
      return (new SingularValueDecomposition(this).norm2());
   }

   /** Infinity norm
   @return    maximum row sum.
   */

   public double normInf ()
   {
      double f = 0;
      for (int i = 0; i < m; i++)
      {
         double s = 0;
         for (int j = 0; j < n; j++)
         {
            s += Math.Abs(A[i][j]);
         }
         f = Math.Max(f,s);
      }
      return f;
   }

   /** Frobenius norm
   @return    sqrt of sum of squares of all elements.
   */

   public double normF ()
   {
      double f = 0;
      for (int i = 0; i < m; i++) 
      {
         for (int j = 0; j < n; j++)
         {
             f = Maths.hypot(f,A[i][j]);
         }
      }
      return f;
   }

   /**  Unary minus
   @return    -A
   */

   public Matrix uminus () 
   {
      Matrix X = new Matrix(m,n);
      double[][] C = X.getArray();
      for (int i = 0; i < m; i++) 
      {
         for (int j = 0; j < n; j++)
         {
            C[i][j] = -A[i][j];
         }
      }
      return X;
   }

   /** C = A + B
   @param B    another matrix
   @return     A + B
   */

   public Matrix plus (Matrix B)
   {
      checkMatrixDimensions(B);
      Matrix X = new Matrix(m,n);
      double[][] C = X.getArray();
      for (int i = 0; i < m; i++) 
      {
         for (int j = 0; j < n; j++)
         {
            C[i][j] = A[i][j] + B.A[i][j];
         }
      }
      return X;
   }

   /** A = A + B
   @param B    another matrix
   @return     A + B
   */

   public Matrix plusEquals (Matrix B) {
      checkMatrixDimensions(B);
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            A[i][j] = A[i][j] + B.A[i][j];
         }
      }
      return this;
   }

   /** C = A - B
   @param B    another matrix
   @return     A - B
   */

   public Matrix minus (Matrix B) {
      checkMatrixDimensions(B);
      Matrix X = new Matrix(m,n);
      double[][] C = X.getArray();
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            C[i][j] = A[i][j] - B.A[i][j];
         }
      }
      return X;
   }

   /** A = A - B
   @param B    another matrix
   @return     A - B
   */

   public Matrix minusEquals (Matrix B) {
      checkMatrixDimensions(B);
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            A[i][j] = A[i][j] - B.A[i][j];
         }
      }
      return this;
   }

   /** Element-by-element multiplication, C = A.*B
   @param B    another matrix
   @return     A.*B
   */

   public Matrix arrayTimes (Matrix B) {
      checkMatrixDimensions(B);
      Matrix X = new Matrix(m,n);
      double[][] C = X.getArray();
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            C[i][j] = A[i][j] * B.A[i][j];
         }
      }
      return X;
   }

   /** Element-by-element multiplication in place, A = A.*B
   @param B    another matrix
   @return     A.*B
   */

   public Matrix arrayTimesEquals (Matrix B) {
      checkMatrixDimensions(B);
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            A[i][j] = A[i][j] * B.A[i][j];
         }
      }
      return this;
   }

   /** Element-by-element right division, C = A./B
   @param B    another matrix
   @return     A./B
   */

   public Matrix arrayRightDivide (Matrix B) {
      checkMatrixDimensions(B);
      Matrix X = new Matrix(m,n);
      double[][] C = X.getArray();
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            C[i][j] = A[i][j] / B.A[i][j];
         }
      }
      return X;
   }

   /** Element-by-element right division in place, A = A./B
   @param B    another matrix
   @return     A./B
   */

   public Matrix arrayRightDivideEquals (Matrix B) {
      checkMatrixDimensions(B);
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            A[i][j] = A[i][j] / B.A[i][j];
         }
      }
      return this;
   }

   /** Element-by-element left division, C = A.\B
   @param B    another matrix
   @return     A.\B
   */

   public Matrix arrayLeftDivide (Matrix B) {
      checkMatrixDimensions(B);
      Matrix X = new Matrix(m,n);
      double[][] C = X.getArray();
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            C[i][j] = B.A[i][j] / A[i][j];
         }
      }
      return X;
   }

   /** Element-by-element left division in place, A = A.\B
   @param B    another matrix
   @return     A.\B
   */

   public Matrix arrayLeftDivideEquals (Matrix B) {
      checkMatrixDimensions(B);
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            A[i][j] = B.A[i][j] / A[i][j];
         }
      }
      return this;
   }

   /** Multiply a matrix by a scalar, C = s*A
   @param s    scalar
   @return     s*A
   */

   public Matrix times (double s) {
      Matrix X = new Matrix(m,n);
      double[][] C = X.getArray();
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            C[i][j] = s*A[i][j];
         }
      }
      return X;
   }

   /** Multiply a matrix by a scalar in place, A = s*A
   @param s    scalar
   @return     replace A by s*A
   */

   public Matrix timesEquals (double s) {
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            A[i][j] = s*A[i][j];
         }
      }
      return this;
   }

   /** Linear algebraic matrix multiplication, A * B
   @param B    another matrix
   @return     Matrix product, A * B
   @exception  IllegalArgumentException Matrix inner dimensions must agree.
   */

   public Matrix times (Matrix B) {
      if (B.m != n) {
         throw new IllegalArgumentException("Matrix inner dimensions must agree.");
      }
      Matrix X = new Matrix(m,B.n);
      double[][] C = X.getArray();
      double[] Bcolj = new double[n];
      for (int j = 0; j < B.n; j++) {
         for (int k = 0; k < n; k++) {
            Bcolj[k] = B.A[k][j];
         }
         for (int i = 0; i < m; i++) {
            double[] Arowi = A[i];
            double s = 0;
            for (int k = 0; k < n; k++) {
               s += Arowi[k]*Bcolj[k];
            }
            C[i][j] = s;
         }
      }
      return X;
   }

   /** LU Decomposition
   @return     LUDecomposition
   @see LUDecomposition
   */

   public LUDecomposition lu () {
      return new LUDecomposition(this);
   }

   /** QR Decomposition
   @return     QRDecomposition
   @see QRDecomposition
   */

   public QRDecomposition qr () {
      return new QRDecomposition(this);
   }

   /** Cholesky Decomposition
   @return     CholeskyDecomposition
   @see CholeskyDecomposition
   */

   public CholeskyDecomposition chol () {
      return new CholeskyDecomposition(this);
   }

   /** Singular Value Decomposition
   @return     SingularValueDecomposition
   @see SingularValueDecomposition
   */

   public SingularValueDecomposition svd () {
      return new SingularValueDecomposition(this);
   }

   /** Eigenvalue Decomposition
   @return     EigenvalueDecomposition
   @see EigenvalueDecomposition
   */

   public EigenvalueDecomposition eig () {
      return new EigenvalueDecomposition(this);
   }

   /** Solve A*X = B
   @param B    right hand side
   @return     solution if A is square, least squares solution otherwise
   */

   public Matrix solve (Matrix B) {
      return (m == n ? (new LUDecomposition(this)).solve(B) :
                       (new QRDecomposition(this)).solve(B));
   }

   /** Solve X*A = B, which is also A'*X' = B'
   @param B    right hand side
   @return     solution if A is square, least squares solution otherwise.
   */

   public Matrix solveTranspose (Matrix B) {
      return transpose().solve(B.transpose());
   }

   /** Matrix inverse or pseudoinverse
   @return     inverse(A) if A is square, pseudoinverse otherwise.
   */

   public Matrix inverse () {
      return solve(identity(m,m));
   }

    /// <summary>
    /// Moore Penrose PseudoInverse
    /// </summary>
    /// <returns></returns>
    public Matrix PInverse(Matrix h)
    {
        SingularValueDecomposition svd = h.svd();

        Matrix V = svd.getV();
        Matrix S = svd.getS();
        Matrix U = svd.getU();

        Matrix U_T = U.transpose();

        Matrix S_T = S.transpose();
        Matrix S_TS = S_T.times(S);
        Matrix S_TS_1 = S_TS.inverse();

        Matrix inversa = V.times(S_TS_1).times(S_T).times(U_T);
        return inversa;
    }

   /** Matrix determinant
   @return     determinant
   */

   public double det () {
      return new LUDecomposition(this).det();
   }

   /** Matrix rank
   @return     effective numerical rank, obtained from SVD.
   */

   public int rank () {
      return new SingularValueDecomposition(this).rank();
   }

   /** Matrix condition (2 norm)
   @return     ratio of largest to smallest singular value.
   */

   public double cond () {
      return new SingularValueDecomposition(this).cond();
   }

   /** Matrix trace.
   @return     sum of the diagonal elements.
   */

   public double trace () {
      double t = 0;
      for (int i = 0; i < Math.Min(m,n); i++) {
         t += A[i][i];
      }
      return t;
   }

   /** Generate matrix with random elements
   @param m    Number of rows.
   @param n    Number of colums.
   @return     An m-by-n matrix with uniformly distributed random elements.
   */

   public static Matrix random (int m, int n) {
      Random random = new Random();
      Matrix A = new Matrix(m,n);
      double[][] X = A.getArray();
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            X[i][j] = random.NextDouble();
         }
      }
      return A;
   }

   /** Generate identity matrix
   @param m    Number of rows.
   @param n    Number of colums.
   @return     An m-by-n matrix with ones on the diagonal and zeros elsewhere.
   */

   public static Matrix identity (int m, int n) {
      Matrix A = new Matrix(m,n);
      double[][] X = A.getArray();
      for (int i = 0; i < m; i++) {
         for (int j = 0; j < n; j++) {
            X[i][j] = (i == j ? 1.0 : 0.0);
         }
      }
      return A;
   }


   /** Print the matrix to stdout.   Line the elements up in columns
     * with a Fortran-like 'Fw.d' style format.
   @param w    Column width.
   @param d    Number of digits after the decimal.
   */

   //public void print (int w, int d) {
   //   print(new PrintWriter(System.out,true),w,d); }

   /** Print the matrix to the output stream.   Line the elements up in
     * columns with a Fortran-like 'Fw.d' style format.
   @param output Output stream.
   @param w      Column width.
   @param d      Number of digits after the decimal.
   */

   //public void print (PrintWriter output, int w, int d) {
   //   DecimalFormat format = new DecimalFormat();
   //   format.setDecimalFormatSymbols(new DecimalFormatSymbols(Locale.US));
   //   format.setMinimumIntegerDigits(1);
   //   format.setMaximumFractionDigits(d);
   //   format.setMinimumFractionDigits(d);
   //   format.setGroupingUsed(false);
   //   print(output,format,w+2);
   //}

   /** Print the matrix to stdout.  Line the elements up in columns.
     * Use the format object, and right justify within columns of width
     * characters.
     * Note that is the matrix is to be read back in, you probably will want
     * to use a NumberFormat that is set to US Locale.
   @param format A  Formatting object for individual elements.
   @param width     Field width for each column.
   @see java.text.DecimalFormat#setDecimalFormatSymbols
   */

   //public void print (NumberFormat format, int width) {
   //   print(new PrintWriter(System.out,true),format,width); }

   // DecimalFormat is a little disappointing coming from Fortran or C's printf.
   // Since it doesn't pad on the left, the elements will come out different
   // widths.  Consequently, we'll pass the desired column width in as an
   // argument and do the extra padding ourselves.

   /** Print the matrix to the output stream.  Line the elements up in columns.
     * Use the format object, and right justify within columns of width
     * characters.
     * Note that is the matrix is to be read back in, you probably will want
     * to use a NumberFormat that is set to US Locale.
   @param output the output stream.
   @param format A formatting object to format the matrix elements 
   @param width  Column width.
   @see java.text.DecimalFormat#setDecimalFormatSymbols
   */

   //public void print (PrintWriter output, NumberFormat format, int width) {
   //   output.println();  // start on new line.
   //   for (int i = 0; i < m; i++) {
   //      for (int j = 0; j < n; j++) {
   //         String s = format.format(A[i][j]); // format the number
   //         int padding = Math.Max(1,width-s.Length()); // At _least_ 1 space
   //         for (int k = 0; k < padding; k++)
   //            output.print(' ');
   //         output.print(s);
   //      }
   //      output.println();
   //   }
   //   output.println();   // end with blank line.
   //}

   /** Read a matrix from a stream.  The format is the same the print method,
     * so printed matrices can be read back in (provided they were printed using
     * US Locale).  Elements are separated by
     * whitespace, all the elements for each row appear on a single line,
     * the last row is followed by a blank line.
   @param input the input stream.
   */

   //public static Matrix read (BufferedReader input) /*throws java.io.IOException*/ {
   //   StreamTokenizer tokenizer= new StreamTokenizer(input);

   //   // Although StreamTokenizer will parse numbers, it doesn't recognize
   //   // scientific notation (E or D); however, Double.valueOf does.
   //   // The strategy here is to disable StreamTokenizer's number parsing.
   //   // We'll only get whitespace delimited words, EOL's and EOF's.
   //   // These words should all be numbers, for Double.valueOf to parse.

   //   tokenizer.resetSyntax();
   //   tokenizer.wordChars(0,255);
   //   tokenizer.whitespaceChars(0, ' ');
   //   tokenizer.eolIsSignificant(true);
   //   java.util.Vector v = new java.util.Vector();

   //   // Ignore initial empty lines
   //   while (tokenizer.nextToken() == StreamTokenizer.TT_EOL);
   //   if (tokenizer.ttype == StreamTokenizer.TT_EOF)
   // throw new java.io.IOException("Unexpected EOF on matrix read.");
   //   do {
   //      v.addElement(Double.valueOf(tokenizer.sval)); // Read & store 1st row.
   //   } while (tokenizer.nextToken() == StreamTokenizer.TT_WORD);

   //   int n = v.size();  // Now we've got the number of columns!
   //   double row[] = new double[n];
   //   for (int j=0; j<n; j++)  // extract the elements of the 1st row.
   //      row[j]=((Double)v.elementAt(j)).doubleValue();
   //   v.removeAllElements();
   //   v.addElement(row);  // Start storing rows instead of columns.
   //   while (tokenizer.nextToken() == StreamTokenizer.TT_WORD) {
   //      // While non-empty lines
   //      v.addElement(row = new double[n]);
   //      int j = 0;
   //      do {
   //         if (j >= n) throw new java.io.IOException
   //            ("Row " + v.size() + " is too long.");
   //         row[j++] = Double.valueOf(tokenizer.sval).doubleValue();
   //      } while (tokenizer.nextToken() == StreamTokenizer.TT_WORD);
   //      if (j < n) throw new java.io.IOException
   //         ("Row " + v.size() + " is too short.");
   //   }
   //   int m = v.size();  // Now we've got the number of rows.
   //   double[][] A = new double[m][];
   //   v.copyInto(A);  // copy the rows out of the vector
   //   return new Matrix(A);
   //}
#endregion

   /** Check if size(A) == size(B) **/

   private void checkMatrixDimensions (Matrix B) 
   {
      if (B.m != m || B.n != n)
      {
         throw new IllegalArgumentException("Matrix dimensions must agree.");
      }
   }


    }
}
